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1.0  INTRODUCTION 

This  memorandum  addresses  issues  involving  gas-phase  finite  rate  chemistry  and  droplet 
vaporization  as  simulated  by  the  LESLIE3D  computer  program.  LESLIE3D  is  an  acronym  for 
Large  Eddy  Simulation  with  Linear  Eddy  model  in  3  Dimension  and  is  developed  by  Professor 
Suresh  Menon  at  the  Georgia  Institute  of  Technology.  Although  this  computer  program  is  the 
product  of  ongoing  research,  it  is  adept  at  simulating  the  multiphase  physics  of  mixtures  of  gases 
with  particles  or  droplets.  Finite  rate  chemical  reaction  mechanisms  are  incorporated  within  this 
computer  program  for  an  arbitrary  number  of  species  and  reactions  (constrained  only  by  available 
computer  memory).  LESLIE3D  has  powerful  capabilities  for  the  simulation  of  gas  phase 
chemistry.  Surface  chemistry  for  solid  particles  is  also  an  inherent  capability  of  the  code,  e.g.,  the 
combustion  of  Aluminum  particles.  Simulation  of  the  evaporation  and  subsequent  heating  and 
combustion  of  liquid  droplets  is  a  historical,  organic  capability  for  this  computer  program.  These 
algorithms  spring  from  one  of  LESLIE3D’s  primary  development  efforts,  the  simulation  of 
turbulent  flow  and  chemistry  in  gas  turbine  combustors. 

1.1  Core  LESLIE3D  Algorithms 

No  discussion  of  LESLIE3D  is  complete  without  a  brief  exposition  of  LESLIE3D’s  core 
numerical  algorithms.  Although  there  are  many  algorithms  undergirding  the  code’s  representation 
of  continuum  dynamics,  three  major  algorithms  stand  out  from  the  others.  First,  LESLIE3D  has 
state-of-the-art  capabilities  for  the  Large  Eddy  Simulation  (LES)  of  compressible  turbulence. 
Suresh  Menon’s  Locally  Dynamic  subgrid  Kinetic  energy  Model  (LDKM)  is  directly  incorporated 
in  the  code.[l]  To  review,  the  idea  behind  LES  is  simple.  Organized  fluid  motions  can  be  ordered 
in  terms  of  scale,  large  size  eddies  cascading  down  to  small  size  eddies.  Informally,  LES 
establishes  a  line  of  demarcation  between  the  large  and  small  scales. [2]  Eddies  within  the  large 
range  of  scales  are  resolved  numerically  by  solving  the  filtered  Navier-Stokes  equations. [3]  The 
properties  of  small  scale  eddies  are  determined  by  modeling.  The  assumption  implicitly  made  is 
that  at  the  small  scales,  turbulent  motions  are  universal.  [4]  Various  mathematical  models  are  used 
to  compute  the  small  scale  (subgrid)  properties  based  upon  properties  of  the  main  flow  in  time.  As 
its  principal  driver,  LDKM  exploits  the  similarity  between  the  Leonard  stress  tensor  at  the  so- 
called  test  scale  with  the  subgrid  scale  stress  tensor.  On  the  axis  of  turbulent  scales,  the  test  scale 
is  a  finite  locus  of  scales  that  are  both  numerically  resolved  and  modeled.  Without  getting  into  the 
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mathematical  details,  similarity  between  these  two  scales  requires  use  of  the  subgrid  and  subtest 
scale  kinetic  energy  for  the  turbulence.  From  non-zero  initial  values,  the  behavior  of  subgrid 
kinetic  energy  is  described  by  an  evolution  equation  that  is  included  in  the  system  of  governing 
equations. [4]  This  combination  of  numerical  simulation  and  modeling  culminates  in  the 
production  of  time  and  space-dependent  coefficients  for  the  production  and  dissipation  of 
turbulence  in  space.  Hence,  LDKM  is  truly  a  local  LES  model. 

The  second  major  capability  of  LESLIE3D  is  its  capacity  for  capturing  the  physics  of 
shock-turbulence.  Shock  waves  are  most  often  modeled  as  discontinuities  in  the  flow  field.  That 
is  to  say,  across  the  shock  wave  front,  flow  properties  such  as  pressure,  density,  temperature, 
entropy  and  individual  velocity  components  change  abruptly.  The  Euler  equations,  a  truncated  set 
of  Navier-Stokes  equations  (disregarding  viscous  effects),  are  often  used  to  model  shock  waves 
since  the  Euler  equations  admit  discontinuous  solutions.  Unfortunately,  discontinuous  solutions 
are  not  admitted  by  the  Navier-Stokes  equations.  Even  for  tiny  regions  near  the  shock  front, 
solutions  must  be  smooth.  In  order  to  automatically  resolve  shock  waves,  the  use  of  upwind  or 
asymmetrical  differencing  is  required.  [5]  On  the  other  hand,  to  resolve  turbulence,  centered  or 
symmetrical  differencing  is  required.  This  situation  presents  a  dilemma.  The  asymmetrical  stencils 
required  for  shock  wave  resolution  possess  a  great  deal  of  numerical  viscosity;  the  associated 
dissipation  tends  to  damp  subtle  turbulent  fluctuations  from  the  flow  field.  Conversely,  the 
simulation  of  turbulence  favors  the  use  of  symmetric,  high  order  computational  stencils.  By  design, 
symmetric  (or  centered)  stencils  incorporate  information  from  both  directions  defined  by  the 
stencil.  For  shocked  flow  fields,  the  bidirectional  incorporation  of  information  induces  instability 
in  the  shock  wave  solver.  The  instability  manifests  itself  in  violent  oscillations  of  flow  field 
properties  originating  near  the  shock  locus  and  spreading  outward.  These  oscillations  ultimately 
corrupt  the  numerical  solution,  induce  instability  and  cause  it  to  diverge. 

LESLIE3D  remedies  this  dilemma  in  two  ways.  First,  we  draw  upon  the  mathematics  of 
the  Piece-wise  Parabolic  Method  (PPM)  to  implement  a  very  low  dissipation  shock-capturing 
scheme.  [6]  The  reconstruction  process  mimics  the  “flattening”  procedure  used  by  PPM  to  limit 
the  amount  of  dissipation  used  to  grant  a  smooth  shock  solution.  A  variant  of  the  Harten,  Lax  and 
van  Leer  Contact-preserving  shock-capturing  scheme  is  employed  along  with  the  Einfeld 
correction  method  for  carbuncle  instabilities  (HLLC/E).[7,8]  For  the  preservation  of  turbulent 
fluctuations,  this  scheme  performs  on  par  with  PPM.  The  second  and  somewhat  more  complex 
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algorithm  implemented  for  shock-turbulence  capture  is  referred  to  as  the  “hybrid”  scheme.  The 
hybrid  scheme  utilizes  a  “switch”  to  apply  HLLC/E  in  the  neighborhood  of  shock  waves  and  a 
high  order  centered  MacCormack  scheme  in  smooth  regions  of  the  flow  field.  This  method 
performs  quite  well  and  has  been  validated  in  number  of  studies.  [4] 

The  third  major  capability  for  LESLIE3D  is  provided  by  Lagrangian  particle/droplet 
tracking  algorithms  for  both  diffuse  and  dense  dispersed  phases.  These  particle  methods  are  fully 
coupled  with  the  gas  phase  to  provide  for  the  highly  accurate  simulation  of  fuel  combustor  flow 
fields.  In  its  original  incarnation,  LESLIE3D  could  only  simulate  diffuse  particle  fields  where  the 
particle  volume  fraction  is  rather  small.  In  recent  years,  the  dispersed  field  algorithms  have  been 
upgraded  to  address  dense  collections  of  particles. [9]  The  difference  lies  in  the  fact  that  dense 
particle  fields  exert  an  internal  “pressure”  on  themselves.  Particles  or  droplets  entrained  in  diffuse 
fields  are  not  subjected  to  this  pressure.  For  dense  dispersed  fields,  this  pressure  is  very  important 
in  determining  the  forces  exerted  on  particles.  This  inter-particle  force  is  not  found  in  diffuse  fields. 
Instead,  diffuse  fields  do  include  the  drag  force.  This  force  is  determined  through  the  use  of 
empirical  drag  laws.  [10] 

The  information  presented  above  provides  only  the  briefest  introduction  to  some  of 
LESLIE3D’s  capabilities  from  the  viewpoint  of  physics,  not  computer  programming.  In  the 
discussions  that  follow,  we  focus  on  LESLIE3D’s  turbulent,  finite  rate  chemistry  algorithms. 
These  algorithms  accurately  capture  gas  phase  reaction  kinetics,  but  more  importantly,  the  physics 
of  turbulent  mixing  is  also  accurately  represented.  It  is  very  important  to  understand  that  for 
turbulent  flow  fields,  if  the  turbulent  motions  are  not  accurately  simulated,  the  rates  for  chemical 
reactions  are  incorrectly  computed,  and  inaccurate  numerical  solutions  result. 

1.2  Combustion  of  Vapor  Clouds 

A  class  of  scenarios  generating  widespread  engineering  interest  involves  the  combustion 
of  a  cloud  of  vapor  released  into  the  atmosphere.  [  1 1 , 1 2]  For  our  purposes,  the  constituents  of  cloud 
may  be  thought  of  in  the  classical  sense  of  combustion.  That  is  to  say,  the  cloud  is  composed  of 
the  fuel  (a  combustible  chemical,  something  that  likes  to  give  up  electrons)  and  the  oxidizer  (a 
substance  that  likes  to  accept  electrons).  Combustible  fuels  always  present  safety  concerns.  These 
fuels  are  usually  stored  in  tanks  (either  above  or  below  ground),  and  given  that  many  of  these  tanks 
are  quite  large  and  can  store  a  great  deal  of  chemical  energy.  It  must  also  be  realized  that  this  type 
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of  fuel  only  bums  when  in  the  gas  phase.  Before  a  liquid  fuel  can  bum,  it  must  first  vaporize.  Note 
that  a  partially  filled  liquid  fuel  tank  contains  an  ullage  volume  (a  void  space  above  the  surface  of 
the  liquid  fuel  volume).  If  this  tank  is  heated  (even  by  the  sun),  part  of  the  fuel  may  vaporize  into 
the  ullage  volume  and  subsequently  pressurize  the  tank.  A  relief  valve  or  vent  installed  on  the  tank 
can  vent  fuel  vapor  into  the  atmosphere  and  present  a  combustion  hazard.  If  no  such  valve  is 
present,  then  it  is  possible  that  the  over  pressurized  tank  may  rupture  releasing  large  amounts  of 
fuel  into  the  surroundings.  This  scenario  constitutes  a  significant  combustion  hazard.  If  ignited  by 
some  random  circumstance,  then  the  fuel  may  begin  to  bum  atomizing  and  vaporizing  more  fuel. 
If  unchecked,  this  combustion  process  can  cascade  into  a  deflagration  or  detonation.  Note  also  that 
light  hydrocarbon  fuels,  e.g.,  acetylene  and  propane  exist  as  gases  at  room  temperature,  so  they 
are  ready  to  combust  if  initiated.  For  these  reasons,  the  combustion  of  fuel  vapor  clouds  deserves 
study  in  the  interest  of  engineering  safer  fuel  storage  and  handling  systems.  Illustrating  some 
aspects  of  the  study  is  the  focus  of  this  memorandum  with  a  particular  emphasis  on  turbulent 
chemistry. 

2.0  Technical  Aspects  of  Turbulent,  Gas  Phase  Combustion 

This  memorandum  is  fundamental  in  that  it  focuses  on  gas  phase  chemistry  in  a  turbulent 
(and  in  some  cases)  shocked  flow  environment.  Gas  phase  chemistry  is  exclusively  considered  in 
order  to  establish  the  physical  framework  for  chemical  reactions.  The  involvement  of  fuel  droplets 
along  with  fuel  atomization  and  evaporation  is  also  very  important,  but  these  multiphase  processes 
inevitably  lead  to  the  fundamental  process  of  gas  phase  combustion. 

2.1  Finite  Rate  Chemical  Mechanisms 

A  finite  rate  chemical  mechanism  consists  of  a  set  of  distinct  chemical  reactions  involving 
a  set  of  chemical  species.  The  term  “finite  rate”  is  used  because  each  chemical  reaction  requires  a 
finite,  non-zero  amount  of  time  to  run  to  completion.  For  this  simple  example,  the  chemical  species 
are  identified  as  compounds  or  elements.  Consider  the  combustion  of  acetylene  gas.  Acetylene 
(C2H2)  is  a  simple  hydrocarbon  fuel  characterized  as  an  alkyne  since  the  Carbon  atoms  are  joined 
by  a  triple  covalent  bond.  A  diagram  of  this  molecule  is  shown  in  Figure  1 . 
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H-C=C-H 

Figure  1.  Molecular  Structure  of  Acetylene 

Acetylene  is  quite  reactive  due  to  the  presence  of  the  triple  bond  since  this  bond  implies  that  three 
pairs  of  electrons  are  being  shared  between  the  Carbon  atoms.  Carbon  has  four  electrons  in  its 
outer  shell.  Since  it  would  like  to  complete  this  shell,  it  would  like  to  share  four  other  electrons.  It 
does  so,  in  this  case,  by  sharing  three  electrons  with  its  neighboring  Carbon,  and  one  electron  with 
the  nearest  Hydrogen  atom.  The  electron  cloud  formed  by  the  shared  electrons  of  the  triple  bond 
is  an  irresistible  temptation  for  Oxygen,  an  element  that  loves  to  acquire  electrons.  The 
fundamental  combustion  reaction  for  acetylene  can  be  written  as  follows. 

C2H2  +  5/2  02  -»•  2  C02  +  H20  (2.1.1) 

For  those  who  may  not  be  familiar  with  chemical  reaction  equations,  the  species  listed  on  the  left 
side  of  the  equation  are  called  reactants;  the  species  occurring  on  the  right  side  of  the  equation  are 
called  products.  Complete  combustion,  as  implied  by  Equation  (2.1.1),  does  not  occur  in  reality. 
In  truth,  other  chemical  products  are  formed  in  small  quantities,  especially  when  a  disproportionate 
amount  of  fuel  is  available  compared  with  oxidizer.  Still,  acetylene  is  a  very  light  hydrocarbon 
(molecular  weight  26.02  g/mol),  so  with  a  preponderance  of  oxygen,  we  expect  a  high  level  of 
conversion  into  the  products  specified  by  (2.1.1).  The  coefficients  (1,  5/2,  2  and  1)  shown  for  the 
species  in  (2.1.1)  are  termed  as  stoichiometric  coefficients.  Although  these  coefficients  are  easily 
interpreted  as  dimensionless  quantities,  they  in  fact  have  the  dimension  of  molecules  (a  cardinal 
number)  or  in  groups  as  moles.  Recall  that  a  mole  contains  Avogardro's  number  of  (6.023  x  1023) 
molecules.  In  this  light,  one  mole  of  acetylene  can  combust  with  5/2  moles  of  oxygen  gas  to 
produce  two  moles  of  Carbon  dioxide  and  one  mole  of  water  vapor.  Of  course,  realistic  problems 
involve  scalar  multiples  of  (2.1.1). 

Equation  (2.1.1)  represents  a  static  concept;  it  confers  no  information  on  the  rate  at  which 
the  reaction  progresses.  Rate  information  is  provided,  in  one  view,  by  the  conventional  theory  of 
chemical  kinetics.  In  a  more  general  case,  reactions  such  as  (2.1.1)  have  both  forward  (r/)  and 
backward  (n)  rates. [13]  These  rates  may  be  formulated  as  follows: 

rf=kf\  C2H2]*1[02r  (2.1.2) 
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n  =  kb  [C02f 1  [IhOf2 


(2.1.3) 


Parameters  kf  and  kb  are  the  forward  and  backward  rate  coefficients.  The  square  bracket  notation 
[  ]  indicates  a  measure  of  concentration,  generally  moles  per  cubic  centimeter  (mol/cm3).  The 
parameters  R\,  Ri,  P\  and  Pi  are  exponents  for  the  concentration  factors.  Under  the  conventional 
theory,  these  parameters  along  with  components  of  the  rate  coefficients  must  be  determined  from 
laboratory  experiments.  More  specifically,  special  curve  fitting  routines  are  required.  Often,  as  a 
first  approximation,  these  exponents  are  chosen  as  the  matching  stoichiometric  coefficients.  But 
for  more  detailed  calibrations  of  the  kinetics,  widely  different  values  are  chose  for  these  exponents. 
The  general  form  for  a  rate  coefficient  may  be  written  as 

k=ATnex  p(-jjj)  (2.1.4) 

Equation  (2.1.4)  is  cast  in  Arrhenius  form,  an  approximation  of  the  potential  energy  surface  for 
the  reaction.  Specifically,  A  is  denoted  as  the  pre-exponential  factor;  T  is  the  absolute  temperature, 
and  Ea  is  the  activation  energy. [13]  The  activation  energy  is  the  “height”  of  the  potential  energy 
barrier  that  must  be  overcome  for  the  reaction  to  proceed.  On  inspection,  (2.1.4)  is  intuitive;  note 
that  for  higher  values  of  the  activation  energy,  the  reaction  rate  is  reduced.  Higher  temperatures 
are  required  to  compensate  for  the  increase  in  Ea.  Also,  for  exothermic  (energy  liberating)  reactions 
such  as  acetylene  combustion,  the  reaction  rate  is  proportional  to  temperature,  so  an  increase  in 
temperature  elevates  the  reaction  rate.  The  final  term  in  (2. 1 .4),  R,  is  the  universal  gas  constant. 

As  it  happens,  for  acetylene  combustion,  the  backward  rate  is  quite  small,  so  we  consider 
only  the  forward  reaction  (2.1.1)  with  the  rate  (2.1.2).  The  data  required  by  (2.1.4)  is  provided  in 
Table  1. 


Table  1.  Arrhenius  Rate  Coefficient  Data  for  Single-Step  Acetylene  Combustion 


A  (mol/cm3/s) 

Ea  (kcal/mol) 

Xi 

*2 

6.5  x  1012 

30 

0.5 

1.25 

This  acetylene  reaction  kinetics  model  and  data  is  used  for  an  example  problem  presented  later  in 
this  memorandum. 
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2.2  The  Closure  Problem  for  Turbulent  Chemistry 

The  Locally  Dynamic  subgrid  Kinetic  energy  Model  (LDKM)  evolves  from  the  filtered 
Navier-Stokes  equations.  [4]  As  a  result,  new  terms  appear  in  the  equations;  these  terms  involve 
averages  of  moments  of  flow  properties.  This  mathematical  process  differs  from  the  classical 
Reynolds  Averaging  method  to  some  degree  as  it  establishes  the  requisite  separation  in  scales  of 
motion.  These  new  terms  have  the  effect  of  adding  variables  to  the  system  of  equations  so  that  the 
total  number  of  variables  exceeds  the  total  number  of  equations.  At  face  value,  this  system  is 
unsolvable;  hence,  the  closure  problem  results.  It  is  necessary  to  add  equations  to  the  system  and 
achieve  mathematical  closure.  For  the  purposes  of  this  work,  we  consider  only  the  species 
evolution  equation,  i.e., 


^  +  +  Y,7  +  0<!f]  =  A  (2-2-1) 

In  this  equation,  p  is  the  local  density  for  the  species  mixture;  Yu  is  the  mass  fraction  of  species  k, 
and  the  Vijc  are  diffusion  velocity  components  for  the  kh  species.  Specifically,  the  filtered  version 
of  this  expression  is 


V.  - 

i,k  Yk  dXi 


(2.2.2) 


where  Dk  is  the  diffusion  coefficient  for  species  k.  The  terms  added  by  the  filtering  process  are 
Ytf ,  0-f  and  in  another  sense,  d)k,  the  filtered  reaction  rate.  The  first  two  terms  are  the  subgrid 


convective  species  flux  and  the  subgrid  diffusive  species  flux  6-  f  .  The  formal  expressions 
for  these  terms  are 

Y,7  =  Pe((u,Yk)  -  u,fk )  (2.2.3) 


eif  =  p{k,A)  ~  Vl, A) 


(2.2.4) 


In  (2.2.3)  and  (2.2.4),  the  angle  brackets  indicate  the  same  mass-average  denoted  by  the  ~  symbol. 
In  practice,  the  subgrid  convective  species  flux  is  modeled;  on  the  other  hand,  the  subgrid  diffusive 
species  flux  is  neglected  as  it  is  thought  to  be  small  compared  to  other  terms.  The  more  critical 
term  for  turbulent  chemistry  is  the  filtered  reaction  rate.  For  this  term,  the  Eddy  Break-Up  (EBU) 
model,  a  simple  turbulent  closure,  is  applied.  [14] 
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The  EBU  model  compares  two  separate  reaction  rates,  first,  the  kinetic  rate  produced  by 
the  chemical  mechanism  and  second,  a  mixing  rate  that  takes  turbulent  motions  into  account.  For 


the  latter  rate,  the  mixing  time  scale  is  defined  as 

_  CEBU-at 

Tmix  “  V2W 


(2.2.5) 


In  this  equation,  CEBU  is  constant,  usually  set  at  one.  The  solution  time  step  is  At,  and  ksgs  is  the 
space  and  time  dependent  subgrid  kinetic  energy  determined  by  the  evolution  equation 


=  £-[(pvt  +  p) 


—  1  +  a 


pd 


(M, 


sgs Y  (pSfesgs) 

f  '  l  Dksss  ) 


d/es§s 

+ 

pvtR  df 

dxi 

Prt  dxi 

'sgs 

duj 

dxt 

+  pcE 

;1 


(2.2.6) 


Equation  (2.2.6)  is  integrated  along  with  the  governing  Navier-Stokes  equations.  The  turbulent 
mixing  rate  is  inversely  proportional  to  this  mixing  time  by  a  factor  related  to  the  driving 
stoichiometric  coefficients.  The  EBU  reaction  rate  is  chosen  as  the  minimum  of  the  kinetic  and 
turbulent  mixing  reaction  rates.  In  most  cases,  the  magnitude  of  the  kinetic  rate  far  exceeds  the 
turbulent  rate.  If  the  kinetic  rate  is  applied  in  turbulent  flow  simulations,  the  chemical  reaction 
rates  are  over  predicted  by  a  significant  margin.  The  classical  EBU  model  differs  from 
LESLIE3D’s  model,  and  the  differences  are  worthy  of  a  brief  description.  The  main  assumption 
behind  EBU  is  that  chemistry  does  not  play  an  explicit  role;  instead,  turbulent  motions  govern  the 
reaction  rate.  In  this  realization,  the  gas  phase  reaction  zone  is  comprised  of  volumes  of  both 
unbumed  and  burned  gases  mixed  by  turbulent  eddies.  The  average  reaction  rate  is  directly 
proportional  to  the  level  of  temperature  fluctuations  and  inversely  proportional  to  a  turbulence 
time  scale.  The  level  of  temperature  fluctuations  is  given,  in  part,  by  the  expression 

0  =  (2.2.7) 

t2-t±  v  7 

as  is  based  upon  laminar  flame  theory  [15]  where 

T2  =  T1+^f  (2.2.8) 

In  these  formulas,  Q  is  the  chemical  energy  liberated  by  the  combustion  reaction;  Yf  is  the  mass 
fraction  of  fuel,  and  Cp  is  the  mixture  constant  pressure  specific  heat.  On  the  other  hand,  the 
classical  EBU  time  scale  is  modeled  as 
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k 


(2.2.9) 


tEBU  — 


£ 


where  k  and  s  are  the  kinetic  energy  and  energy  dissipation  rate  for  the  turbulence.  With  the  use  of 


these  quantities,  the  classical  EBU  reaction  rate  is  defined  as 


<*>6  ~  Cebu p-0(l  —  O') 


(2.2.10) 


In  this  equation,  p  is  the  mixture  density  and  the  tilde  symbol  indicates  an  average.  Of  course, 
Cebu  is  a  constant  chosen  for  the  model.  LESLIE3D’s  EBU  model  is  implemented  differently.  An 
EBU  reaction  rate  is  computed  for  each  reaction.  All  of  the  species  used  in  the  problem  (reactants 
and  products)  are  concatenated  in  order.  If  a  particular  species  is  a  product  (non-zero 
stoichiometric  coefficient  I )■),  then  set 


A  special  caveat  on  (2.2. 1 1)  is  that  if  the  kinetic  rate  for  the  reaction  goes  forward,  then  it  depends 


solely  on  reactant  concentrations,  so  we  set  P  equal  zero.  Conversely,  if  the  species  is  a  reactant 
(non-zero  stoichiometric  coefficient  Rj ),  we  set 


(2.2.12) 


If  the  reaction  proceeds  in  the  reverse  direction,  it  depends  solely  on  the  product  concentrations, 
so  R  is  set  to  zero  in  this  case.  With  the  use  of  these  details,  the  LESLIE3D  EBU  reaction  rate  is 
computed  as 


P—R 


(2.2.13) 


for  the  reaction  in  question.  The  EBU  model,  in  either  of  the  forms  shown  here,  is  an  effective  and 
relatively  inexpensive  closure  for  turbulent  chemistry.  Other  methods  certainly  exist  since  closures 
are,  in  fact,  models. 

2.3  Initializing  the  Turbulent  Field 

Although  Menon’s  LDKM  methodology  constitutes  a  state-of-the-art  model,  the  subgrid 
kinetic  energy  equation  is  an  evolution  equation  that  requires  a  non-zero  initial  turbulent  velocity 
and  subgrid  kinetic  energy  field  in  order  to  simulate  the  evolution  of  the  turbulence.  For  any  three- 
dimensional  volume,  there  are  an  infinite  number  of  possible  turbulent  velocity  distributions;  we 
must  select  one  such  distribution  for  initializing  LESLIE3D.  For  problems  involving  shock  waves, 
we  may  define  a  non-zero  velocity  and  subgrid  kinetic  energy  field  at  the  shock  wave  outflow. 
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Since  the  resolved  velocity  field  is  non-zero,  the  turbulence  cannot  decay  before  the  simulation 
begins.  This  idea  is  rather  ad  hoc  since  we  do  not  know  a  priori  the  turbulence  velocity  field  in 
this  region.  For  problems  involving  initially  stagnant  (zero  velocity),  shock-free  flow  fields,  the 
aforementioned  procedure  does  not  apply.  Instead,  one  must  define  a  field  of  isotropic  turbulence 
consistent  with  the  governing  equations  in  the  ambient  gas  volume.  For  even  non-shocked  flow 
fields,  this  turbulence  can  evolve  for  many  solution  iterations  without  decaying,  and  it  is 
completely  physical.  Moreover,  this  field  is  defined  strictly  in  terms  of  three-dimensional,  resolved 
field  velocity  components;  the  subgrid  kinetic  energy  is  initialized  consistent  with  the  resulting 
flow  field  kinetic  energy.  Below,  the  mathematical  methodology  is  described  for  initializing  the 
initial  turbulent  velocity  field.  [16] 

Begin  by  considering  a  nearly  quiescent  flow  field  with  velocity  field  v  =  v(x1,x2,x3). 
For  the  flow  field  to  be  incompressible,  we  require  that  V  ■  v  =  0.  By  performing  a  Fourier  integral 
transform  on  this  equation,  one  obtains 

v-k  =  0  (2.3.1) 

where  the  wavevector  is  defined  as 

k  =  (kx,ky,kz)  (2.3.2) 

The  transformed  variables  are  indicated  by  the  symbol.  With  (2.3.1)  in  mind,  suppose  that  we 
start  with  a  more  general  velocity  field  v  where  (2.3.1)  does  not  hold  true.  It  is  necessary  to  alter 
this  velocity  field  so  that  (2.3.1)  does  hold.  For  the  general  v,  its  multi-dimensional  Fourier 
transform  may  be  written  as  follows. 

v  =  %  +  (2.3.3) 

Since  (2.3.3)  is  tediously  annotated,  what  this  equation  really  says  is  that  a  the  Fourier 
transformation  of  a  general  velocity  field  can  be  decomposed  into  the  sum  of  velocity  fields 

parallel  and  perpendicular  to  wavevector  k.  From  what  this  equation  says,  in  the  light  of  (2.3.1), 
the  velocity  field  that  we  need  is  the  one  perpendicular  to  k,  so  we  solve  (2.3.3)  for  this  part  of  the 
transformed  velocity  field,  i.e., 

=  $  ~  V\\%  (2-3.4) 

Observing  the  sense  of  the  vector  k,  we  note  that  is  that  component  of  v  in  the  direction  of 
the  unit  vector  k  where 
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Hence, 


k 

k 


(2.3.5) 


(2.3.6) 


By  substituting  into  (2.3.4),  we  obtain 


(2.3.7) 


For  computational  purposes,  it  is  useful  to  rewrite  (2.3.7)  in  indicial  notation  where 


(2.3.8) 


Equation  (2.3.8)  is  the  formula  needed  to  create  an  incompressible  velocity  field  from  an  arbitrary 
non-conservative  velocity  field. 

To  compute  the  incompressible  velocity  field,  as  mentioned  above,  we  start  with  an 
arbitrary  velocity  distribution  produced  by  using  a  pseudorandom  number  generator  to  create  the 
velocity  components.  [17]  In  this  case,  the  magnitude  of  each  random  velocity  vector  is  kept  fairly 
small.  The  next  step  is  to  perform  a  discrete  Fourier  transform  on  this  velocity  field.  Each  block  in 
the  computational  flow  field  is  fully  structured;  that  is,  each  point  in  the  grid  is  indexed  by  the 
ordered  triple  (/',  /,  k)  where  1  <  i  <  imax,  1  <  j  <  jmax  and  1  <  k  <  kmax.  To  support  the  loop 
structure  of  the  transform,  it  is  helpful  to  require  that  imax,  jmax  and  kmax  be  odd  natural  numbers. 
Then,  the  field  cell  counts  in  the  index  directions  become 


icmax  =  imax  -  1 
jcmax  =  jmax  -  1 
kcmax  =  kmax  -  1 


(2.3.9) 


The  counting  index  maxima  for  the  transform  become 


Ni  =  icmax/2 
Ni  =  jcmax/2 
7V3  =  kcmax/2 


(2.3.10) 


Let  the  velocity  vector  be  denoted  as 


V  =  (u-L.^.U-z) 


(2.3.11) 


then  the  forward  transform  may  be  written  as 
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where  l  —  iVp/2  <  kp  <  Np/2,  p=l,2,3.  As  a  matter  of  review,  (2.3. 12)  relies  upon  the  Euler 
formula 

exp(i0)  =  cos(0)  +  i  sin(0)  (2.3.13) 

Note  that  the  i  in  (2.3.12)  is  the  imaginary  unit.  Moreover,  the  computed  Fourier  coefficients  dpare 
complex  numbers.  That  is  to  say,  while  the  velocity  components  up  are  numbers  on  the  real  line, 
the  Fourier  coefficients  tip  exist  in  the  complex  plane.  These  coefficients  are  ready  for  the 
incompressibility  correction. 

The  correction  of  the  up  to  enforce  incompressibility  (more  pertinently,  mass 
conservation),  is  implemented  in  the  same  three  nested  loops  implied  by  the  transform  (2.2.12). 
Instead  of  this  transform,  equation  (2.3.8)  is  calculated.  In  pseudo-code,  the  loops  may  be 
expressed  as  follows. 

Fori  -  NJ2  <kt<  NJ 2,1  -  N2/ 2  <k2<  N2/ 2,1  -  N3/2  <  k3  <  JV3/ 2,  with  kp± 

0,  for  all  p  =  1,2,3 ,  we  compute 

Re(f  •  k )  =  Re(u1)k1  +  Re(u2)k2  +  R e(u3)k3  (2.3.14) 

Im(f  •  k)  =  Im^)/^  +  hn(u2)k2  +  Im (u3)k3  (2.3.15) 

k2  =  kl  +  kl  +  kl  (2.3.16) 

The  real  and  imaginary  parts  of  the  transformed  velocity  components  are  corrected  by  the 
statements  below. 

Re(Ui)  =  Re(Ui)  -  ^Re(d  •  k)  (2.3.17) 

Im^)  =  Im^)  —  ^lm(v  ■  k)  (2.3.18) 

Re(u2)  =  Re(u2)  -  ^Re(d  •  k)  (2.3.19) 

Im(u2)  =  Im(u2)  —  ^|lm(v  •  k )  (2.3.20) 

Re(u3)  =  Re(u3)  -  Re(^  •  k)  (2.3.21) 

Im(u3)  =  Im(u3)  —  ^|lm(v  •  k )  (2.3.22) 

It  is  important  to  realize  that  the  transformed  velocity  components  must  be  set  equal  zero  if  all 
three  wavenumber  components  are  zero.  If  this  action  is  not  taken,  a  spurious,  non-physical  mode 
is  left  within  the  transform.  This  undesirable  mode  is  analogous  to  a  rigid  body  mode  occurring  in 
finite  element  structural  analysis. 
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The  transformed  velocity  components  now  represent  an  incompressible  flow  field 
satisfying  mass  conservation.  However,  this  field  is  not  representative  of  any  particular  turbulent 
spectrum.  That  is  to  say,  isotropic  turbulence  has  a  particular  kinetic  energy  spectrum  (when 
computed  versus  the  scalar  wavenumber  locus).  Since  the  wavenumber  field  is  defined  in  terms  of 
a  discrete  subset  of  the  natural  numbers,  then  it  has  a  well-defined  maximum  and  a  finite  number 
of  elements.  Here,  the  wavevector  is  addressed  by  its  modulus  squared,  a  scalar  locus.  By  using 
(2.3.10),  the  maximum  squared  wavenumber  is 

fcLx  =  Ni2  +  IVf  +  N$  (2.3.23) 

Not  all  of  the  natural  numbers  between  zero  and  /c„iax  represent  acceptable  wavevectors.  For  these 
values  of  k2,  the  kinetic  energy  is  assigned  the  default  value  of  zero.  For  valid  values  of  k2  in 
(2.3.16),  the  kinetic  energy  is  computed  as  follows. 

£(k2)  =  i  (lk2=£.£  ( u(k )  u(k)  +  v{k)v{k)  +  w*(k)w(/c)))  (2.3.24) 


The  raw  spectrum  is  obtained  by  summing  the  energy  over  the  different  wavenumber  bands.  The 
different  wavenumber  bands  can  be  determined  by  dividing  k2nax  by  max(Ni,N2,Ns).  In  each  band, 
the  spectral  magnitude  is  scaled  by  a  factor  to  achieve  the  magnitude  of  an  archived  spectrum  such 
as  Comte-Bellot  and  Corrsin.[18]  This  scale  factor  is  directly  applied  to  the  transformed  velocity 
coefficient. 

The  final  step  in  developing  an  isotropic  field  is  to  apply  the  inverse  discrete  Fourier 
transformation  to  the  coefficients  derived  above.  The  inverse  transformation  is  written  as 


u(n1,n2,n3)  =  77^77  Zf1  1 


Ni N2N3  ^k2=i  J^k3=i_^3 


u  (/cx,  /c2,  A:3) 


■  exp| i^k-tri-L  +  k2n2  +  k3n3)]  (2.3.25) 

forO  <  rii  <  /V;  —  1,  Uj  =  1,2,3.  For  uniform,  Cartesian  geometries,  this  velocity  distribution 
is  representative  of  an  isotropic  turbulence  field. 


3.0  TEST  PROBLEMS 

Due  to  the  myriad  of  physics  involved  in  turbulent  chemistry,  phase  transitions  and 
atomization  for  fuel  droplets,  simulating  the  combustion  process  for  liquid  fuel  is  highly  complex. 
This  complexity  is  inherent  in  both  computer  code  development  and  problem  execution  and 
solution.  In  an  attempt  to  this  describe  type  of  simulation  in  a  practical  way,  a  set  of  test  problems 
is  introduced  below.  In  this  work,  the  calorically  perfect  gas  equation  of  state  is  utilized  for  the 
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test  problems. 


3.1  Gas  Phase  Acetylene  Combustion 

This  problem  is  designed  to  illustrate  gas  phase  combustion;  no  liquid  droplets  are 
involved.  Still,  the  scenario  is  of  practical  interest  for  the  field  of  fuel  storage  safety.  Acetylene  is 
a  highly  reactive,  explosive  fuel  commonly  used  as  a  welding  gas.  Since  it  is  a  gaseous  substance 
at  atmospheric  temperature  and  pressure,  it  must  be  stored  in  pressurized  tanks.  Should  this  type 
of  tank  rupture,  the  pressurized  contents  may  ignite.  A  simplified  version  of  this  problem  can  be 
simulated  by  defining  an  air-acetylene  distribution  near  a  solid  surface.  The  distribution  of  species 
is  designed  to  vary  from  a  premixed  gas  volume  at  the  exterior  of  the  distribution  to  the  diffusion 
gas  volume  near  the  distribution  core.  The  species  involved  in  this  problem  are:  C2H2  (acetylene), 
N2  (nitrogen),  O2  (oxygen),  Ar  (argon),  CO2  (carbon  dioxide)  and  H2O  (water).  The  single  step 
reaction  mechanism  presented  in  equation  (2.1.1)  is  applied  here.  Naturally,  nitrogen  and  argon 
are  non-reactive  species,  but  since  they  have  noticeable  mass  fractions  in  air,  they  are  included  to 
make  the  mixing  problem  more  realistic.  The  mass  fractions  for  the  species  in  dry  air  are  presented 
in  Table  2. 


Table  2.  Mass  Fractions  of  Chemical  Species  in  Air 


n2 

02 

Ar 

C02 

0.7552 

0.2314 

0.0129 

0.0005 

The  computational  domain  is  box-shaped  10  meters  long  by  10  meters  wide.  In  terms  of  height, 
the  vertical  domain  is  stretched  in  the  region  above  5  meters  to  accommodate  an  expanding  cloud. 
The  computational  grid  is  comprised  of  48  blocks;  each  block  has  20  cells  in  each  of  three 
directions.  The  fuel  cloud  is  centered  at  x  =  z  =  5  meters;  y  =  0.  The  acetylene,  nitrogen  and  oxygen 
mass  fractions  are  described  by  a  Gaussian  distribution  applied  to  the  acetylene  mass  fraction.  The 
relative  mass  fractions  of  nitrogen  and  oxygen  match  those  in  air.  The  mass  fraction  of  acetylene 
is  given  by 

W*)  =  (3.1.1) 
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Figure  2.  Iso-surfaces  of  the  equivalence  ratio  for  the  initial  acetylene  distribution  (oblique  view) 


Figure  3.  Iso-surfaces  of  the  equivalence  ratio  for  the  initial  acetylene  distribution  (overhead  view) 
where  xcen  contains  the  fuel  cloud  center  coordinates,  and  x  is  a  point  in  the  domain.  For  this  case, 
the  width  of  the  Gaussian  is  implicitly  defined  by  the  parameter  a.  For  this  study,  a  is  set  equal 
one.  An  equivalence  ratio  distribution  can  be  computed  for  the  acetylene-air  cloud.  For  the 
initial  conditions,  iso-surfaces  of  the  equivalence  ratio  are  shown  in  Figure  2.  An  oblique  view  is 
shown  here.  A  symmetric  view  from  above  is  shown  in  Figure  3.  The  combustion  reaction  is 
initiated  by  placing  a  high  temperature  distribution  within  the  computational  domain.  This 
distribution  is  also  Gaussian  in  the  form  of  (3.1.1)  with  o  =  0.25.  The  center  of  the  distribution  is 
shifted  to  the  coordinates  x  =  z  =  3.75;  y  =  0.  This  configuration  is  easily  simulated  with 
LESLIE3D;  the  LES  and  EBU  algorithms  are  activated  for  turbulent  combustion.  The  EBU  case 
is  mixing  dominated  with  comparably  slow  reaction  progress.  Iso-surfaces  of  equivalence  ratio  are 
shown  for  this  problem  at  0.83  seconds  in  Figure  4  and  compared  with  the  initial  equivalence 
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Figure  5.  Iso-surfaces  of  equivalence  ratio  for  the  acetylene  distribution  at  6.36  seconds  (EBU  model) 


Figure  6.  Iso-surfaces  of  equivalence  ratio  for  the  acetylene  distribution  at  9.53  seconds  (EBU  model) 


ratio.  It  is  evident  that  the  gas  volume  is  moving  outward  radially,  and  the  fuel  density  is  lowering 
as  is  evidenced  by  the  falling  equivalence  ratio.  Peculiarities  in  the  turbulent  flow  field  cause  the 
bum  front  to  evolve  with  an  odd  shape.  This  effect  is  shown  in  Figure  5  at  6.36  seconds  and  in 
Figure  6  at  9.53  seconds.  As  is  exhibited  in  Figures  7  through  9,  the  evolution  of  the  temperature 
field  is  more  interesting.  Turbulent  motions  are  evident  through  the  motion  of 
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Figure  7.  Iso-surfaces  of  temperature  for  acetylene  combustion  at  0.83  seconds  (EBU  model) 
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Figure  8.  Iso-surfaces  of  temperature  for  acetylene  combustion  at  2.01  seconds  (EBU  model) 


Figure  9.  Iso-surfaces  of  temperature  for  acetylene  combustion  at  9.53  seconds  (EBU  model) 


the  flame  front  in  this  sequence  of  plots.  Although  the  bum  is  initiated  at  single  point,  it  does  not 
spread  in  a  uniform  radial  directory.  Perturbations  in  the  flow  field  caused  by  the  turbulence  split 
the  flame  into  parts  that  propagate  in  different  directions.  Although  the  bum  is  more  intense  near 
the  point  of  initiation,  it  also  develops  two  distinct  lobes  when  viewed  from  above  and  takes  on 
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Figure  10.  Iso-surface  plot  of  subgrid  kinetic  energy  at  9.53  seconds 


Figure  11.  Iso-surface  plot  of  water  mass  fraction  at  9.53  seconds 
the  shape  of  a  bow.  The  flame  spreads  laterally  but  unevenly  for  this  reason  as  is  evident  in  Figures 
7,  8  and  9.  The  spatial  extent  of  the  flame  correlates  well  with  the  distribution  of  subgrid  kinetic 
energy.  One  may  see  this  correlation  by  comparing  Figures  9  and  10.  Figure  10  contains  an  iso¬ 
surface  plot  of  subgrid  kinetic  energy  viewed  perpendicular  to  the  y  axis.  Due  to  the  lack  of  grid 
refinement  used  in  the  case,  the  turbulent  wall  interactions  are  over  predicted.  Still,  the  correlation 
between  the  combustion  reaction  and  subgrid  kinetic  energy  is  well  illustrated.  A  view  of  product 
formation  is  also  provided  by  the  numerical  solution  and  shown  in  Figure  11.  The  concentration 
of  water  closely  follows  the  temperature  and  subgrid  kinetic  energy 
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Figure  12.  Iso-surface  plots  of  temperature  at  0.05  seconds  (kinetic  model) 


Figure  13.  Iso-surface  plots  of  subgrid  kinetic  energy  at  0.07  seconds  (kinetic  model) 
distributions.  The  reaction  is  driven  by  the  chemical  rate  (represented  by  temperature)  and  by 
turbulent  mixing  (represented  by  subgrid  kinetic  energy). 

It  is  interesting  to  contrast  the  results  of  the  EBU  model  with  those  associated  with  is 
referred  to  here  as  the  kinetic  model.  For  the  kinetic  model,  the  reaction  rates  produced  by  the 
chemical  mechanism  are  chosen  as  the  actual  reaction  rates.  No  turbulent  closure  is  applied  to  the 
chemical  rates.  As  a  result,  reactants  found  within  a  grid  cell  are  assumed  to  be  fully  mixed  and 
react  to  the  fullest  extent  dictated  by  the  kinetic  mechanism.  Figures  12  and  13  shown  the  results 
for  temperature  at  0.05  seconds  and  for  subgrid  kinetic  energy  at  0.07  seconds,  respectively,  for 
this  model.  Turbulence  still  has  an  effect  on  the  overall  distribution  of  properties,  but  it  is  less 
pronounced  since  it  has  less  of  an  effect  on  local  reaction  rates.  The  kinetic  model  can  over  predict 
reaction  rates  by  neglect  turbulent  mixing.  Under  the  kinetic  model,  reactions  proceed  more 
rapidly.  This  effect  is  illustrated  by  Figure  14,  a  comparison  of  the  temperature  distribution  for  the 
kinetic  model  at  0.05  and  0.07  seconds.  The  different  temperature  distributions  can  be  identified 
by  the  blue  lines.  It  is  evident  that  the  blast  cloud 
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Figure  14.  Temperature  contours  for  the  kinetic  model  at  0.05  and  0.07  seconds 


Figure  15.  Temperature  contours  for  the  EBU  model  at  0.04  and  0.07  seconds 
expands  significantly  in  the  span  of  0.02  seconds.  A  similar  plot  is  provided  for  the  EBU  model  in 
Figure  15.  The  temperature  contours  closely  plot  over  one  another  in  this  with  0.03  seconds 
between  the  temperature  solutions.  The  kinetic  solution  advances  faster  than  the  EBU  solution. 

3.2  Kerosene  Droplet  Combustion 

The  acetylene  bum  problem  discussed  in  Section  3.1  is  fundamental  in  that  it  entails  the 
combustion  of  a  purely  gaseous  substance.  Acetylene  is  a  purely  gaseous  substance  in  this  case, 
so  no  phase  change  is  involved.  A  more  common  occurrence  entails  the  release  of  fuel  in  droplet 
form;  under  the  right  thermal  conditions,  the  droplets  break  down  and  vaporize.  The  portion  of  the 
fuel  now  in  the  gas  phase  can  combust  transforming  into  reaction  products.  The  test  fuel  considered 
here  is  kerosene,  a  principal  component  of  jet  fuel.  This  problem  is  more  interesting  and  complex 
than  the  preceding  test  case  because  kerosene  is  not  a  pure  substance;  rather,  it  is  mixture  of  alkanes 
(like  heptane),  naphthenes  (or  cyclo-alkanes,  like  cyclobutane),  alkenes  (like  ethylene)  and 
aromatic  hydrocarbons  (like  benzene). [19]  Secondly,  the  evaporation  phase  change  and 
subsequent  vapor  ignition  is  caused  by  the  impact  of  a  series  of  incident  and  reflected  shock  waves 
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inside  of  a  square  shock  tube.  [20]  Heating  of  the  droplets,  in  this  situation,  requires  hydrodynamic 
compression  of  the  gas  with  the  increase  in  temperature  predicted  by  the  equation  of  state.  The 
repetition  of  archived  experiments  via  simulation  is  also  complicated  by  the  fact  that  in  many 
experiments,  the  initial  conditions  in  the  shock  tube,  (e.g.,  pressures  in  the  driver  and  driven 
sections)  are  not  documented.  Experimentalists  often  concentrate  on  achieving  a  particular 
pressure  behind  the  reflected  shock  wave. [21]  The  task  of  replicating  this  pressure  in  the 
simulation  set-up  is  left  to  the  computational  physicist. 

Given  the  time  allowed  for  accomplishing  this  report,  a  simplified  problem  involving  the 
shock  tube  ignition  of  kerosene  droplets  is  described.  The  shock  tube  has  a  square  cross-section 
that  is  one  meter  in  width  by  one  meter  in  height.  The  tube  length  is  48  meters,  and  the  driver 
section  is  16  meters  long  (1/3  of  the  tube  length).  To  generate  a  reflected  shock  wave,  the  driven 
end  of  the  tube  is  terminated  with  a  solid  surface.  Each  one  meter  section  of  tube  length  is 
designated  as  a  subdomain  (or  block).  Each  block  is  uniformly  meshed  with  20  cells  in  each  of  the 
three  Cartesian  directions.  The  field  of  20,000  kerosene  droplet  parcels  is  suspended  between 
locations  24  and  33.6  meters  along  the  length  of  the  shock  tube.  Each  parcel  consists  of  100 
droplets.  The  particle  radii  are  randomly  selected  to  be  between  37.5  pm  and  50.0  pm.  The  droplet 
(x,  y,  z)  coordinates  are  given  by  the  equations 

xd  =  24.0  +  9.6  •  RAND 

yd  =  0.2  +  0.6  •  RAND  (3.2.1) 

zd  =  0.2  +  0.6  •  RAND 

These  equations  place  the  droplets  at  random  locations  within  a  box-shaped  region  within  the  tube. 
RAND  is  a  series  of  sequentially  generated  pseudorandom  numbers  used  to  place  the  droplets  at 
random  locations.  By  using  this  process,  we  avoid  creating  a  Cartesian  array  of  droplets.  Also,  the 
droplets  are  assigned  small,  random  initial  velocities  (magnitude  0.1  m/s)  since  most  experiments 
inject  droplets  into  the  shock  tube  with  some  initial  velocity. 

LESLIE3D  contains  a  model  for  evaporation  that  draws  upon  a  thermochemical  properties 
file  for  kerosene.  As  a  result,  this  computer  program  can  simulate  the  evaporation  of  kerosene 
droplets.  As  the  droplets  begin  to  evaporate  into  kerosene  vapor,  finite  rate  chemistry  algorithms 
are  invoked  to  bum  the  kerosene  into  products.  A  simplified  version  of  the  reduced  reaction 
kinetics  model  developed  by  Franzelli  et  al.  [22]  is  employed  to  capture  the  chemistry  in 
conjunction  with  the  EBU  turbulent  closure  model.  This  mechanism,  consisting  of  two  steps,  is 
very  efficient  for  numerical  applications.  It  is  written  as 
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(3.2.2) 


(1)  C10H20  +  10  02  -►  10  CO  +  10  H20 

The  second  step  involves  the  oxidation  of  carbon  monoxide  to  for  carbon  dioxide,  i.e., 

(2)  CO  +  Vi  O2  C02  (3.2.3) 


Table  3.  Kerosene  Combustion  Reaction  Coefficient  Data 


Reaction  Step 

A 

n 

E/R 

1 

8x10“ 

1 

4.15  x  104 

2 

4.5  xlO10 

1 

2  x  104 

The  rate  expressions  for  reactions  (3.2.2)  and  (3.2.3)  are  written  as  follows. 

rn  =A:1[C10H20]0-55[O2]0-9  (3.2.4) 

rf2=k2[C0f[02f5  (3.2.5) 

Correction  terms  based  upon  equivalence  ratio  are  applied  to  the  pre-exponential  factors  in  [19]; 
for  simplicity,  these  factors  are  not  employed  here.  For  the  version  of  LESLIE3D  addressed  here, 
the  user  is  required  to  write  the  source  code  for  any  particular  finite  rate  chemistry  mechanism 
unless  such  coding  is  available  from  an  archived  source. 

Initial  conditions  within  the  driver  and  driven  sections  of  the  shock  tube  are  provided  in 
Table  4.  A  turbulent  initial  velocity  field  is  also  established  within  the  shock  tube  in  accordance 
with  the  method  described  in  Section  2.3.  The  conditions  are  selected  to  ensure  droplet  evaporation 
and  ignition  via  shock  heating.  The  initial  data  for  both  the  gas  and  particle  fields  is 


Table  4.  Shock  Tube  Initial  Conditions 


Property 

Driver  Section 

Driven  Section 

Pressure  (Pa) 

214967.0 

2293.14 

Temperature  (K) 

800.0 

300.0 

coded  a  priori  in  initial  LESLIE3D  restart  and  droplet  files.  These  files  are  read  by  LESLIE3D  at 
time  of  code  execution.  By  far,  the  majority  of  time  in  conducting  this  simulation  is  in  setting  up 
the  initial  data;  very  little  time  is  required  to  configure  LESLIE3D.  On  48  cores,  this  simulation 
executes  rather  quickly,  at  approximately  three  seconds  per  iteration  (or  step). 

Results  have  been  generated  for  the  kerosene  droplet  bum  problem  out  to  1 1 1  milliseconds 
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(ms)  over  7000  iterations.  A  sequence  of  temperature  field  plots  is  included.  Each  plot  has  an 
overlay  of  the  kerosene  droplet  field  existing  at  the  same  time.  Figure  16  contains  these  plots  in 
order  at  12  ms,  15  ms,  18  ms,  21  ms,  23  ms  and  25  ms.  In  this  figure,  note  that  although  the 
temperature  field  is  shown  in  translucent  coloration,  the  particle  parcels  are  shown  as  block  dots. 
During  this  segment  of  the  numerical  solution,  kerosene  droplet  parcels  are  actively  evaporating. 
LESLIE3D  actively  provides  the  parcel  count;  a  sampling  of  this  count  is  shown  in  Table  4.  With 
evaporation  combined  with  hydrodynamic  and  turbulent  motion,  the 

Table  3.  Time  dependent  inventory  of  kerosene  droplet  parcels 


Time  (ms) 

12 

15 

18 

21 

23 

25 

Parcels 

20,000 

18,375 

12,641 

6,419 

1,097 

83 

gaseous  kerosene  field  evolves.  Slice  plots  of  the  mass  fraction  for  kerosene  are  presented  in  Figure 
17  at  12  ms,  15  ms  ,  18  ms,  21  ms,  23  ms  and  26  ms.  From  this  figure  one  may  see  that  kerosene 
vapor  is  formed  at  the  shock  outflow  and  is  then  advected  along  with  flow  field  trailing  the  main 
shock.  By  plotting  kerosene’s  mass  fraction  on  the  same  scale,  the  accumulation  of  kerosene 
(increasing  mass  fraction)  is  evident  at  the  shock  outflow.  Since  the  droplets  are  suspended  near 
the  tube  core,  away  from  the  walls,  little  kerosene  vapor  accumulates  near  the  tube  walls.  For  this 
reason,  the  slice  plots  are  taken  along  a  plane  through  the  tube  centerline.  Doing  so  provides  a 
clearer  picture  of  the  kerosene  vapor  distribution.  This  emission  of  kerosene  vapor  motivates  the 
question  as  to  whether  it  is  burning  in  the  tube.  An  answer  to  this  inquiry  is  provided  by  examining 
the  generation  of  product  mass,  e.g.,  for  water  vapor,  within  the  tube.  Recall  that  for  each  mole  of 
kerosene  burned,  the  first  reaction  step  generates  10  moles  of  water  and  10  moles  of  carbon 
dioxide.  Slice  plots  of  the  mass  fraction  of  water  are  shown  in  Figure  1 8.  It  is  interesting  to  compare 
Figures  17  and  18.  The  distribution  of  water  vapor  closely  follows  that  of  kerosene  vapor  because 
kerosene  is  burning  producing  the  water  vapor.  Both  of  these  gases  are  advected  by  the  shock 
outflow;  hence,  as  is  shown  in  these  two  figures,  the  gas  volumes  follow  the  shock  wave.  The 
distribution  of  carbon  monoxide  would  look  very  similar  since  it  is  produced  on  a  one  to  one  molar 
ratio  with  water.  Similar  statements  can  be  made  for  the  distribution  of  carbon  dioxide  in  the  tube. 
The  mass  fraction  is  about  an  order  of  magnitude  smaller  than  that  of  carbon  monoxide  its  primary 
reactant. 
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Figure  16.  Shock  tube  temperature  plots  including  particle  positions  at  12, 15, 18,  21,  23  and  25  milliseconds 
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Figure  17.  Slice  plots  of  kerosene  vapor  mass  fraction  at  12, 15, 18,  21,  23  and  26  milliseconds 


H20 

— -2.41 1  e-06 

( 

- 1 .6079e-6 


1 .2053 e -6 


-6.0264e-7 

-0.000e+00 

H20 

--2.4 11  e-06 

- 1 .6079e-6 


■  I 


1 .20538-6 


-6.0264e-7 

-0.0008400 

H20 

--2.4 11  e-06 

- 1 .6079e-6 


■  1  1.2053e-6 


-6.0264e-7 

* 

■=0.0006400 

H2Q 

--2.4 11  e-06 

-1.60798-6 


M] 


1 .20538-6 


-6.02648-7 

-0.0006+00 


26 

Distribution  A 


H20 

--2.4116-06 

I 

- 1 .6079e-6 


■ 


1 .20536-6 


-6.02646-7 

CO.OOOe+OO 

Figure  18.  Slice  plots  of  water  vapor  mass  fraction  at  15, 18,  21,  23  and  26  milliseconds 


4.0  CONCLUSIONS 

This  technical  memorandum  documents  continued  progress  in  assimilating  the  LESLIE3D 
multiphase  physics  computer  program  for  applications  in  shocked  reactive  flow  fields.  As  an 
expository  work,  this  memorandum  concentrates  on  describing  LESLIE3D’s  turbulent  finite  rate 
chemistry  algorithms.  The  term  “turbulent”  is  used  for  two  reasons.  In  the  first  place,  LESLIE3D 
incorporates  state-of-the-art  algorithms  for  the  large  eddy  simulation  of  compressible  turbulence. 
Secondly,  the  code  applies  turbulent  chemical  closures  to  incorporate  the  physical  limitations  of 
reaction  rates  imposed  by  mixing.  The  simplest  example  of  the  turbulent  chemical  closure  is  the 
Eddy  Break-Up  model  described  herein.  The  effect  of  this  closure  is  illustrated  by  the  acetylene 
combustion  problem.  The  evolution  of  the  acetylene  bum  is  compared  with  and  without  this 
closure.  Results  show  that  the  combustion  reaction  proceeds  much  more  quickly  for  the  kinetic 
model.  This  model  does  not  use  a  turbulent  closure;  rather  in  each  cell,  perfect  mixing  is  assumed 
for  the  reactants.  When  the  Eddy  Break-Up  model  is  employed,  the  reaction  proceeds  more  slowly. 
Since  perfectly  mixing  rarely  occurs  in  the  real  world,  we  expect  the  kinetic  model  to  overestimate 
chemical  reaction  rates.  This  assertion  has  important  ramifications  for  a  number  of  different 
engineering  fields. 

A  second  topic  addressed  by  this  memorandum  is  the  evaporation  and  subsequent 
combustion  of  liquid  fuel  droplets.  Kerosene,  a  complex  hydrogen  mixture,  is  explored  from  the 
standpoint  of  shock-initiated  combustion.  LESLIE3D  simulates  the  evaporation  of  a  wide 
distribution  of  droplet  diameters  and  then  applies  a  finite  rate  chemical  mechanism  to  simulate 
kerosene  combustion.  A  related  fact  is  that  LESLIE3D  can  simulate  chemical  mechanisms 
involving  an  arbitrary  number  of  species  and  reaction  steps  known  a  priori  (limited  only  by 
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computer  memory).  In  this  work,  a  simplified  adaptation  of  the  mechanism  developed  by  Franzetti 
et  al.  is  employed.  The  simulation  results  show  the  time  dependent  evaporation  of  kerosene  (as 
modeled)  and  its  transformation  into  products.  The  positions  and  velocities  of  the  droplets  are 
predicted  by  LESLIE3D  in  the  Lagrangian  sense.  The  successful  completion  of  this  problem 
constitutes  a  highly  applicable  framework  for  solving  an  array  of  problems  involving  the 
combustion  of  fuel  or  in  the  wider  view,  the  consumption  of  an  arbitrary  reactant.  Time  permitting, 
more  complex  models  may  be  constructed  and  solved  to  elucidate  important  aspects  of  chemical 
physics. 
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